%Script to perform a simplified turbulence assessment of the AC_FIRE
%simulation presented in Kiefer et al. (2025) [GMD].
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%PRELIMINARIES%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear
b=pwd;b=strsplit(b,'/');b=b{2};
if(strcmp(b,'home'))%102 server (CPU and memory-intensive processing and/or figure generation)
    addpath(genpath('/home/mmfire/matlab_functions/FUNCTIONS'))
    addpath(genpath('/home/mmfire/matlab_functions/SHAPEFILES'))
    arps_dir='/d102c/mmfire/ARPS-CANOPY-DEVS-FIRE-devel/ARPS-CANOPY-DEVS-FIRE_5.4.2_DEVEL_output/ARPS-CANOPY/RUNS/FL1_WS3_UI1/co1_cs1_cp0/lp601_dp3e4/';
    devs_dir='/d102c/mmfire/ARPS-CANOPY-DEVS-FIRE-devel/ARPS-CANOPY-DEVS-FIRE_5.4.2_DEVEL_output/DEVS-FIRE/RUNS/FL1_WS3_UI1/co1_cs1_cp0/lp601_dp3e4/';
    grid_dir='/home/mmfire/ARPS-CANOPY-DEVS-FIRE-devel/ARPS-DEVS_IDEALIZED_TESTS/PHASE2/GRIDINFO_ANALYSIS/';%GRID INFO IS IDENTICAL FOR ARPS and ARPS-CANOPY SIMS.
    plot_dir='/home/mmfire/ARPS-CANOPY-DEVS-FIRE-devel/ARPS-CANOPY-DEVS_CODEDEVEL/ARPS-CANOPY-DEVS-FIRE_5.4.2_DEVEL/diff_plots/RUNS/FL1_WS3_UI1/co1_cs1_cp0/lp601_dp3e4/COUPLING_ANALYSIS/SURF_tightview/';
elseif(strcmp(b,'dska') | strcmp(b,'dskd') | strcmp(b,'dske'))%150 server (CPU and memory-intensive processing and/or figure generation)
    addpath(genpath('/dska/mtkiefer/matlab_functions/FUNCTIONS'))
    addpath(genpath('/dska/mtkiefer/matlab_functions/SHAPEFILES'))
    arps_dir='/dska/mtkiefer/ARPS-CANOPY-DEVS-FIRE_REALCASE_TESTS_output/SERDP_2019_ACDF/ARPS-CANOPY/RUNS/TEST4/';
    devs_dir='/dska/mtkiefer/ARPS-CANOPY-DEVS-FIRE_REALCASE_TESTS_output/SERDP_2019_ACDF/DEVS-FIRE/RUNS/TEST4/';
    devsconfig_dir='/dska/mtkiefer/ARPS-CANOPY-DEVS-FIRE_REALCASE_TESTS/SERDP_2019_ACDF/DEVS-FIRE/RUNS/TEST4/';
    grid_dir='/dska/mtkiefer/ARPS-CANOPY-DEVS-FIRE_REALCASE_TESTS/SERDP_2019_ACDF/ANALYSIS/GRID_INFO/';
    plot_dir='/dska/mtkiefer/ARPS-CANOPY-DEVS-FIRE_REALCASE_TESTS/SERDP_2019_ACDF/ANALYSIS/PLOTS/TEST4/';
    mat_dir='/dska/mtkiefer/ARPS-CANOPY-DEVS-FIRE_REALCASE_TESTS/SERDP_2019_ACDF/ANALYSIS/MATFILES/TEST4/';
elseif(strcmp(b,'Volumes'))%my laptop (simpler processing and/or figure generation)
    addpath(genpath('/Users/mmfire/Documents/MATLAB/matlab_functions/FUNCTIONS'))
    addpath(genpath('/Users/mmfire/Documents/MATLAB/matlab_functions/SHAPEFILES'))
    arps_dir='/Volumes/SAMSUNG/ARPS-CANOPY-DEVS-FIRE-devel_macbookair/ARPS-CANOPY-DEVS-FIRE_REALCASE_TESTS/SERDP_2019_ACDF/ANALYSIS/RUNS/ARPS-CANOPY/TEST4/';
    devs_dir='/Volumes/SAMSUNG/ARPS-CANOPY-DEVS-FIRE-devel_macbookair/ARPS-CANOPY-DEVS-FIRE_REALCASE_TESTS/SERDP_2019_ACDF/ANALYSIS/RUNS/DEVS-FIRE/TEST4/';
    devsconfig_dir='/Volumes/SAMSUNG/ARPS-CANOPY-DEVS-FIRE-devel_macbookair/ARPS-CANOPY-DEVS-FIRE_REALCASE_TESTS/SERDP_2019_ACDF/DEVS-FIRE/RUNS/TEST4/';
    grid_dir='/Volumes/SAMSUNG/ARPS-CANOPY-DEVS-FIRE-devel_macbookair/ARPS-CANOPY-DEVS-FIRE_REALCASE_TESTS/SERDP_2019_ACDF/ANALYSIS/GRID_INFO/';
    plot_dir='/Volumes/SAMSUNG/ARPS-CANOPY-DEVS-FIRE-devel_macbookair/ARPS-CANOPY-DEVS-FIRE_REALCASE_TESTS/SERDP_2019_ACDF/ANALYSIS/PLOTS/TEST4/';
    mat_dir='/Volumes/SAMSUNG/ARPS-CANOPY-DEVS-FIRE-devel_macbookair/ARPS-CANOPY-DEVS-FIRE_REALCASE_TESTS/SERDP_2019_ACDF/ANALYSIS/MATFILES/TEST4/';
end
arps_head='acdf_serdp_0030';


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 6/26/26: READ IN VARIABLES FROM ARPS HISTORY FILES: U, V, W, E(SGS)
% EXTRACT AT AT APPROXIMATE CENTER OF BURN UNIT: I = 124, J = 128
% (39.914902 N 74.594086 W)
% WRITE OUT ALL LEVELS.  WRITE OUT EVERY SECOND.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% jj=128;
% ii=124;
% rst=10620;
% ren=10620+7200;
% rin=1;
% NX=253;NY=253;NZ=83;
% NT=(ren-rst)/rin+1;
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ARPS.U=nan(NT,NZ);ARPS.V=nan(NT,NZ);ARPS.W=nan(NT,NZ);ARPS.E=nan(NT,NZ);
% ARPS.t=nan(NT,1);ARPS.z=nan(1,NZ);%Instantaneous time, unstaggered Zp
% %Loop over files.
% c=1;
% for t=rst:rin:ren
%     disp(num2str(t))
%     if(t<10)
%         numsec=strcat('00000',num2str(t));
%     elseif(t<100)
%         numsec=strcat('0000',num2str(t));
%     elseif(t<1000)
%         numsec=strcat('000',num2str(t));
%     elseif(t<10000)
%         numsec=strcat('00',num2str(t));
%     elseif(t<100000)
%         numsec=strcat('0',num2str(t));
%     else
%         numsec=strcat('',num2str(t));
%     end
%     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%     %mac server:
%     if(c==1)
%         temZP = hdfread(strcat(arps_dir,arps_head,'.hdf',numsec),'/zp','Index',{[1 1 1],[1 1 1],[NZ NY NX]});
%         %compute unstaggered ZP.
%         ZPu=nan(NZ,NY,NX);
%         for k=1:NZ-1
%             ZPu(k,:,:)=0.5*(temZP(k,:,:)+temZP(k+1,:,:));
%         end
%         ZPu(NZ,:,:)=temZP(NZ,:,:);%Same as ARPS unstaggering (in gradsio3d.f90)
%         ARPS.z=double(ZPu(:,jj,ii)'-temZP(2,jj,ii));%HEIGHT ABOVE GROUND LEVEL
%     end
%     temU = hdfread(strcat(arps_dir,arps_head,'.hdf',numsec),'/u','Index',{[1 1 1],[1 1 1],[NZ NY NX]});
%     temV = hdfread(strcat(arps_dir,arps_head,'.hdf',numsec),'/v','Index',{[1 1 1],[1 1 1],[NZ NY NX]});
%     temW = hdfread(strcat(arps_dir,arps_head,'.hdf',numsec),'/w','Index',{[1 1 1],[1 1 1],[NZ NY NX]});
%     temE = hdfread(strcat(arps_dir,arps_head,'.hdf',numsec),'/tke','Index',{[1 1 1],[1 1 1],[NZ NY NX]});
%     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%     %Unstagger the grid (E IS AT CENTER, NO UNSTAGGERING NEEDED)
%     %compute unstaggered U-wind.
%     UWNDu=nan(NZ,NY,NX);
%     for i=1:NX-1
%         UWNDu(:,:,i)=0.5*(temU(:,:,i)+temU(:,:,i+1));
%     end
%     UWNDu(:,:,NX)=temU(:,:,NX);%Same as ARPS unstaggering (in gradsio3d.f90)
%     %compute unstaggered V-wind.
%     VWNDu=nan(NZ,NY,NX);
%     for j=1:NY-1
%         VWNDu(:,j,:)=0.5*(temV(:,j,:)+temV(:,j+1,:));
%     end
%     VWNDu(:,NY,:)=temV(:,NY,:);%Same as ARPS unstaggering (in gradsio3d.f90)
%     %compute unstaggered W-wind
%     WWNDu=nan(NZ,NY,NX);
%     for k=1:NZ-1
%         WWNDu(k,:,:)=0.5*(temW(k,:,:)+temW(k+1,:,:));
%     end
%     WWNDu(NZ,:,:)=temW(NZ,:,:);%Same as ARPS unstaggering (in gradsio3d.f90)
%     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%     %Extract at approximate center of burn unit.
%     ARPS.U(c,:)=UWNDu(:,jj,ii);
%     ARPS.V(c,:)=VWNDu(:,jj,ii);
%     ARPS.W(c,:)=WWNDu(:,jj,ii);
%     ARPS.E(c,:)=temE(:,jj,ii);
%     ARPS.t(c,1)=t;
%     c=c+1;
% end
% clear UWNDu VWNDu WWNDu ZPu temU temV temW temE temZP
% save(strcat(mat_dir,'/forpub/acfire_figures_turbtest.mat'),'ARPS')
% 
% exit

load(strcat(mat_dir,'/forpub/acfire_figures_turbtest.mat'))
ARPS.Es=ARPS.E;ARPS=rmfield(ARPS,'E');%Rename E as Es (SGS)
ARPS.Er=NaN*ARPS.Es;%Initialize resolved TKE array.

%Strategy: Use 1-hour with-fire mean as background to compute perturbations
%(as in Kiefer et al. 2024), then compute instantaneous TKE, then produce
%1-min means (as in Kiefer et al. 2022).  Refer to those mean TKE values in
%the text of the Kiefer et al. (2025) manuscript ("not shown"), and be
%explicit in describing how the numbers were computed (e.g., with-fire
%background, 1-hour means).
NZ=size(ARPS.U,2);NT=size(ARPS.U,1);
for k=1:NZ
    %Compute perturbation U,V,W
    tem = squeeze(ARPS.U(2:NT,k));%"2:": start at 14:57:01 EDT.
    tem_reshaped=reshape(tem,(NT-1)/2,2);
    tem1=cat(1,NaN,reshape(detrend(tem_reshaped,'constant'),NT-1,1));%NaN at 14:57:00 EDT.
    tem = squeeze(ARPS.V(2:NT,k));%"2:": start at 14:57:01 EDT.
    tem_reshaped=reshape(tem,(NT-1)/2,2);
    tem2=cat(1,NaN,reshape(detrend(tem_reshaped,'constant'),NT-1,1));%NaN at 14:57:00 EDT.
    tem = squeeze(ARPS.W(2:NT,k));%"2:": start at 14:57:01 EDT.
    tem_reshaped=reshape(tem,(NT-1)/2,2);
    tem3=cat(1,NaN,reshape(detrend(tem_reshaped,'constant'),NT-1,1));%NaN at 14:57:00 EDT.
    %Compute resolved TKE
    ARPS.Er(:,k)=0.5*((tem1.*tem1)+(tem2.*tem2)+(tem3.*tem3));
    clear tem tem_reshaped tem1 tem2 tem3
end
%Compute 1-min average total TKE.  Summarize distribution.
ARPS.Et=ARPS.Er+ARPS.Es;
Et=nan((NT-1)/60,NZ);Er=nan((NT-1)/60,NZ);Es=nan((NT-1)/60,NZ);

for k=1:NZ
    tem = squeeze(ARPS.Et(2:NT,k));%"2:": start at 14:57:01 EDT.
    tem_reshaped=reshape(tem,60,(NT-1)/60);
    Et(:,k)=mean(tem_reshaped,1);
    tem = squeeze(ARPS.Er(2:NT,k));%"2:": start at 14:57:01 EDT.
    tem_reshaped=reshape(tem,60,(NT-1)/60);
    Er(:,k)=mean(tem_reshaped,1);  
    tem = squeeze(ARPS.Es(2:NT,k));%"2:": start at 14:57:01 EDT.
    tem_reshaped=reshape(tem,60,(NT-1)/60);
    Es(:,k)=mean(tem_reshaped,1);    
end

plot(1:120,squeeze(0.5*(Et(:,4)+Et(:,5))),'ko-',... %6-m AGL/1-min mean/ total TKE
    1:120,squeeze(0.5*(Er(:,4)+Er(:,5))),'ro-',... %6-m AGL/1-min mean/ resolved TKE
    1:120,squeeze(0.5*(Es(:,4)+Es(:,5))),'go-') %6-m AGL/1-min mean/ SGS TKE
legend('Et','Er','Es')
